function elaplace2

%  error solving Laplaces equation as a function of N (M=N and M=N/2)
%  matrix equation solved using CGM 

% clear all previous variables and plots
clear *
clf

% n_points = number of different N values to use
n_points=7;

% loop for various values of N (M=N)
for ic=1:n_points

	% set parameters
	if ic<5
		N=10*2^(ic-1)+2;
	else
		N=100*2*(ic-4)+2;
	end;
	M=N;
	points(ic)=N-2;

	% generate grid
	x=linspace(0,1,N);
	y=linspace(0,1,M);
	h=x(2)-x(1);
	k=y(2)-y(1);
	alpha=h/k;

	% generate matrix and RHS of equation
	A=matrix(alpha,N-2,M-2);
	b=rhs(x,y,alpha,N-2,M-2);

	%tic;
	% use CGM to solve matrix equation
	v=cgm(A,b);
	%toc

	% transform back to u(i,j)
	u=map(x,y,v,N,M);

	% calculate exact solution
	sol=exact(N,M,x,y);

	% calculate error
	error(ic)=max(max(abs(u-sol)));

end;


% loop for various values of N (M=N/2)
for ic=1:n_points

	% set parameters
	if ic<5
		N=10*2^(ic-1)+2;
	else
		N=100*2^(ic-4)+2;
	end;
	M=(N-2)/2+2;
	points2(ic)=N-2;

	% generate grid
	x=linspace(0,1,N);
	y=linspace(0,1,M);
	h=x(2)-x(1);
	k=y(2)-y(1);
	alpha=h/k;

	% generate matrix and RHS of equation
	A=matrix(alpha,N-2,M-2);
	b=rhs(x,y,alpha,N-2,M-2);

	%tic;
	% use CGM to solve matrix equation
	v=cgm(A,b);
	%toc

	% transform back to u(i,j)
	u=map(x,y,v,N,M);

	% calculate exact solution
	sol=exact(N,M,x,y);

	% calculate error
	error2(ic)=max(max(abs(u-sol)));

end;

% get(gcf)
% set(gcf,'Position', [7  477 530 280]);
plotsize(530,280)

% plot solution
loglog(points,error,'-ko','MarkerSize',7)
hold on
grid on
loglog(points2,error2,'-ks','MarkerSize',7)

% define axes and legend used in plot
xlabel('Number of Points Along x-Axis','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
legend(' M = N',' M = N/2',1);

% have MATLAB use certain plot options (all are optional)
set(gca,'MinorGridLineStyle','none')
set(gca,'FontSize',14);
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% exact solution:  assumes gL=gR=gB=0
function sol=exact(N,M,x,y)
sol=zeros(N,M);

if N<100
	modes=200;
else
	modes=400;
end;

for nn=1:modes
	np=nn*pi; e6=exp(6)*(-1)^nn;
	an(nn)=(12/5)*np*( 4*e6*(np^2-36)*(np^2+6) - 672*np^2 - 5*np^4 - 26352)/(36+np^2)^4;
%	an(nn)=an(nn)/sinh(np);
end;

for j=1:M-1
	for i=1:N
		s=0;
		for nn=1:modes
			np=nn*pi; np1=np*(y(j)-1); np2=np*(y(j)+1);
			sinh2=(exp(np1)-exp(-np2))/(1-exp(-2*np));
			s=s+an(nn)*sinh2*sin(np*x(i));
%			s=s+an(nn)*sinh(nn*pi*y(j))*sin(nn*pi*x(i));
		end;
		sol(i,j)=s;
	end;
end;
for i=1:N
	sol(i,M)=gT(x(i));
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% top BC function (y=1)
function g=gT(x)
g=x*(1-x)*(4/5-x)*exp(6*x);

% right BC function (x=1)
function g=gR(y)
g=0;

% bottom BC function (y=0)
function g=gB(x)
g=0;

% left BC function (x=0)
function g=gL(y)
g=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for constructing RHS
function b=rhs(x,y,alpha,N,M)
b=zeros(N*M,1);
for i=1:N
	b(i)=alpha^2*gB(x(i+1))+b(i);
	it=(M-1)*N+i;
	b(it)=alpha^2*gT(x(i+1))+b(it);
end;
for j=1:M
	it=(j-1)*N+1;
	b(it)=gL(y(j+1))+b(it);
	it=j*N;
	b(it)=gR(y(j+1))+b(it);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for creating the matrix A
function A=matrix(alpha,N,M)

nm=N*M;

% generate the diagonal elements
D = sparse(1:nm,1:nm,2*(1+alpha^2)*ones(1,nm),nm,nm);

% generate the subdiagonal elements
SD = sparse(2:nm,1:nm-1,-ones(1,nm-1),nm,nm);
for i=N+1:N:nm-1
	SD(i,i-1)=0;
end;

% generate the sub-subdiagonal elements
SSD=sparse(N+1:nm,1:nm-N,-alpha^2*ones(1,nm-N),nm,nm);

% use symmetry of A to complete construction 
A=D+SD+SD'+SSD+SSD';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for converting back to u(i,j)
function u=map(x,y,v,N,M)
u=zeros(N,M);
for ix=1:N
	u(ix,1)=gB(x(ix));
	u(ix,M)=gT(x(ix));
end;
for iy=1:M
	u(1,iy)=gL(y(iy));
	u(N,iy)=gR(y(iy));
end;
for j=2:M-1
	for i=2:N-1
		l=(j-2)*(N-2)+i-1;
		u(i,j)=v(l);
	end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for the CGM
function x=cgm(A,b)

% set CGM parameters
nm=length(b);
if nm<10000
	tol=0.000001;
else
	tol=0.000000001;
end;
x=zeros(nm,1);

% start iteration
r=b-A*x;
d=r;
rr=dot(r,r);
error=1;
beta=0;
counter=1;
while error>tol
	counter=counter+1;
	q=A*d;
	alpha=rr/dot(d,q);
	x=x+alpha*d;
	r0=r;
	r=r-alpha*q;
	rr0=rr;
	rr=dot(r,r);
	error=norm(alpha*d,inf)/norm(x,inf);
%	error=sqrt(rr);
	beta=rr/rr0;
	d=r+beta*d;
end;

fprintf('\n  Number of CGM Iterations = %i     Error =  %e \n\n',counter,error)

% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);